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1.0 INTRODUCTION 

When working sequences of images, registration of the frames to a common point of reference 
is an essential prerequisite for many types of image analysis. Examples include 

• Coadding images to improve signal-to-noise ratio (SNR) 

• Superresolution 

• Color ratios 

• Polarization ratios 

• Pan-sharpening 

• Conversion of infrared multispectral imagery to emissivity-temperature space 

Example - This section begins with an example showing the importance of accurate frame 
registration. Figure 1 shows an IKONOS satellite image of Venice. The color ratio of an image is 
frequently used for classification in multispectral imagery. The green vs. red color ratio for this 
image was computed by extracting the green and red components of the image and performing 
simple element-by-element matrix division (R1 = Green/Red). The color ratio was computed a 
second time (R2) using a green image with a (deliberate) subpixel registration error of 0.1 pixel. 
Figure 2 displays the quantity (R2-R1)/R1, which is the normalized error in color ratio. Although 
the registration error is small, the color ratio error is as large as 15%. This sensitivity to frame 
registration is typical for many image analysis procedures. Therefore, it is important to choose 
the best frame registration technique available. 




2.0 PHASE CORRELATION METHOD 

The phase correlation method^'^ (PCM) is a popular Fourier domain method to register two 
images. It computes a phase difference map that (ideally) contains a single peak. The location 
of the peak is proportional to the relative translation [dx, dy] between the two images. The PCM 
is resilient to noise and image defects and is readily automated. It is completely equivalent to 
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correlation in the spatial domain, but the calculation is orders of magnitude faster in the Fourier 
domain. The mathematical details are as follows: 

Consider two identical images i1 and i2, with i2 shifted by an amount [Ax, Ay] relative to i1, 

i2(x,y) = il(x-Ax,y-Ay) (1) 

and which obey periodic boundary conditions, 

il{M + x,N + y) = il{x,y) (2) 

where the image size is MxN pixels. Denote Fourier transforms of i1 and i2 by 11 and I2. From 
the Fourier shift theorem, II and 12 differ only by a linear phase term j((jOxAx + WyAy). 
Specifically, 


I2{co^,co^) = n{cD^,co^)e 


y(ffl_,Ax+ a^Ay) 


(3) 


where tOx and cOy are the frequency variables in column and row. The normalized cross-power 
spectrum of the images Cl 2 is defined as 


Cl2(co^,cOy) = 


II. * conj{I2) 
|/1. * conj{I2)\ 


e 


-j{co^Aoc+ cOylAy) 


(4) 


The operator * is the Schur product (also known as the Hadamard element-by-element matrix 
product) and con] is the complex conjugate operator. The right-hand side of Eq. (4) is simply the 
Fourier transform of a Dirac delta function. Stated differently, the inverse Fourier transform of 
the normalized cross-power spectrum Cl 2 is a two-dimensional Dirac delta function 6(x-Ax) 5(y- 
Ay) with a peak location corresponding to the displacement [Ax, Ay] between the two images 
(Fig. 3). Leprince et al.^ show that this result remains unbiased in the presence of additive noise 
or optical blur. The inset in the upper right-hand corner of the phase map shows that for real 
images, which contain both noise and uncorrelated content, the peak is distorted from the ideal 
delta function profile. 
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Figure 3. Shifted Images and IFFT (Cross-Power Spectrum) 
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3.0 ACCURACY TRADE-OFFS 

The accuracy of PCM registration depends upon the content and quality of the images. 
Representative trade-offs of accuracy vs. image noise and array size are indicated in this 
section. Perhaps the most prominent advantage of the PCM is that it is exceptionally robust 
against noise. The accuracy of the PCM for the clean image on the left of Fig. 4 is 0.02 pixels. 
For the severely degraded noisy image on the right of Fig. 4 (11% rms noise), the registration 
accuracy is still 0.15 pixels. Considering the poor quality of the image, this level of accuracy is 
truly remarkable. The accuracy for variable rms noise levels is shown in Fig. 5. The rms noise 
level in Fig. 5 is given in counts, and the image itself has a dynamic range of 8 bits (0 to 255 
counts). Figure 6 shows that, as expected, the registration accuracy improves for larger images. 
The abscissa is the linear dimension of the (square) focal plane array in pixel units. 
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Figure 5. Error vs. Noise 


Figure 6. Error vs. Array Size 


Periodic Boundary Conditions vs. the Reai World - An essential feature of the present work 
is that it uses realistic images that do not obey periodic boundary conditions. The PCM is strictly 
applicable to images that obey periodic boundary conditions (Eq. [2]). However, images of this 
kind are unphysical. Figures 7 and 8 illustrate the unphysical “world view” according to periodic 
boundary conditions. In this “world view,” the original frame repeats endlessly like a mosaic. A 
shifted image (e.g., denoted by the dashed line) contains exactly the same features as the 
original image. The only difference is that the spatial features have been cyclically shifted (Fig. 
9). The result is that the shifted frame can be brought into perfect correlation with the original 
one by a simple cyclic “unshift” of the shifted frame. However, periodic boundary conditions 
hardly 
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ever apply in the real world. Thus a fundamental assumption (Eq. [2]) underlying the PCM is 
hardly ever obeyed by any pair of real images. 


The correct “world view” for two shifted frames is 
illustrated by the yellow and red outlined frames in 
Fig. 9. Each image is a different sample taken from 
the same global scene. The two images share a great 
deal of common features. These are useful for 
correlation, which is at the heart of the PCM. 

However, each frame also contains features near the 
edge that the other does not share. For example, the 
boy is entirely within the image denoted by the red 
frame, but only partially within the image denoted by 
the yellow frame. These uncorrelated features inject 
noise into the PCM correlation. This in turn sets a 
fundamental limit on the accuracy of frame 
registration. If the images are decomposed into correlated vs. uncorrelated parts: 

71 = If + If (5) 



72 = 72^ + 72^ 


( 6 ) 


then Eq. (4) for the normalized cross-product will deviate from the ideal case (delta function) by 
a factor proportional to the fraction of the image that is uncorrelated. 


Error = 


7 ^ 


Area of uncorrelated image 
Area of correlated image 


(7) 


This is an intuitively reasonable result. 

In order to be realistic, the test images used in this work were all generated according to the 
procedure illustrated in Fig. 9. Specifically, all shifted images were selected from a larger global 
scene and contain uncorrelated content near the borders (I.e., as actually happens in the real 
world). They were not generated by cyclic shifting with periodic boundary conditions. This 
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ensures that the performance estimates provided in this work are representative of the 
performance one will actually obtain on real-world images - none of which obey periodic 
boundary conditions. The PCM works in practice because when the shift between the two 
images is small compared to the size of the image, most of the scene is common to both 
frames. Therefore, correlation still works fairly well. In addition, uncorrelated effects at the 
borders of the images are mitigated by multiplying each image by a window function that 
smoothly tapers to zero at the edges of the image. 

Test Images - In addition to nonperiodic boundary conditions, real images are subject to noise 
and frequency aliasing. In addition, some focal plane arrays have dead space between the 
pixels. If one is concerned with registering images to within only ±1 pixel, these real-world 
effects are generally inconsequential. However, they become significant as one attempts to 
push the PCM to subpixel accuracy. Therefore, the test images used to determine algorithm 
performance contain all of these real-world effects. The test images used in this investigation 
(Fig. 10) were derived from high-resolution (approx 3,000 x 3,000 pixels) satellite images of 
terrestrial scenes. The images included diverse scenes - ocean, urban, desert, and rural. 
Downsampled shifted images, such as would be observed by a lower resolution sensor, and 
containing the aforementioned “real-world” effects, were generated as follows. The original high- 
resolution satellite image was first blurred by a Gaussian point-spread function (PSF) to 
simulate the optics blur of the sensor. The blurred image was then downsampled by a factor of 
10 onto a rectilinear focal plane array (FPA) of size NFPA x NFPA, as notionally indicated in 
Fig. 11. The baseline value was NFPA = 256 pixels. The aimpoint of the FPA with respect to the 
(blurred) high-resolution image was offset by integral pixel amounts, as measured in the high 
spatial resolution coordinate system. Thus, after downsampling, the image shifts were exact 


Figure 10. Test Images (Before Blurring) 
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multiples of 0.1 pixel. Each pixel of the simulated FPA 
was surrounded by an optional inactive dead space to 
simulate a fill factor < 100%. (The fill factor of a sensor 
array is the ratio of active area to total area.) The 
nominal fill factor was 0.8'^2 = 64%. Gaussian additive 
noise was added as the last step in the simulation. A 
Gaussian blur kernel was used, so the downsampled 
images have a realistic amount of aliasing. 

Description of Algorithms - Method 1 by Hoge^ uses 
singular value decomposition (SVD) to isolate the 
dominant rank-one principal component of the phase 
difference matrix. This method was developed for 
magnetic resonance imaging. Method 2 by Ren, 

Vlachos, and Jiang"^ (RVJ) uses a subspace projection 
method to fit a linear function (i.e., a plane) to the phase of the normalized cross-power 
spectrum (Eq. [4]). It is the fastest of the three methods. Method 3 by Guizar-Sicairos, Thurman, 
and Fienup^ (GTF) is a traditional PCM approach with a Discrete Fourier Transform interpolation 
routine to accurately determine the peak location with small memory allocation. Each of these 
methods has specific strengths that, if combined into a single code, might possibly outperform 
any of the individual methods. However, this was beyond the scope of this report. Method 3 
turned out to be the most reliable, so two variants of Method 3 were investigated in an effort to 
achieve even better accuracy. Method 4 applies a complex gradient operator to the images® 
before applying Method 3. Specifically, the original image, I, was replaced by the complex 
gradient image, Icg- This operation acts as a high-pass filter to enhance edges and high spatial 
frequencies. 



I 

( 6 ) 

Method 5 applies histogram equalization to the images before applying Method 3. 

Results - This section describes the accuracies of the three registration methods (plus the two 
variants of Method 3) based upon the 11 test images of Fig. 10. The baseline conditions were a 
square focal plane with dimension NFPA = 256 pixels, dynamic range = 8 bit, 10 counts rms 
additive Gaussian noise, and FPA fill factor = 64%. For all methods, the downsampled images 
were apodized by a cosine window function before attempting frame registration. Foorosh et al.^ 
recommend prefiltering of the phase difference matrix to remove aliased components (generally 
at high spatial frequencies). This recommendation was not implemented in the present 
investigation, but it is a good idea that should increase the accuracy of any subpixel PCM 
algorithm. The difficulty is that the filtering must be tailored to each image and sensor. 

Figure 12 and Table 1 show the registration accuracy for Methods 2 through 5. Method 3 (GTF) 
and its two minor variants (Methods 4 and 5) performed best. The performance of Method 2 
(RVJ) was comparable to that of Method 1 (GTF) for 8 out of the 11 images. However, it was 
noticeably poorer for image Nos. 1, 7, and 10. The author could find nothing special about these 
three images to explain why the algorithm of Method 2 performed poorly in these particular 
cases. Method 1 (Hoge) performed poorly for optical imagery, and the results are not shown in 
Fig. 12. Minor improvements to the accuracy of Method 3 were achieved by preprocessing 
according to Methods 4 and 5. Figure 13 shows the robustness vs. additive Gaussian noise 
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(rms counts). Method 2 was the first to “break” (at 50 counts rms noise) under the load of 
increasing noise. This is an extremely high noise level and would not likely be encountered in 
practice. 
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Figure 12. Registration Accuracy 


Figure 13. Noise Sensitivity 


Method 

1 

2 

3 

4 

5 

Mean 

1.210 

0.097 

0.045 

0.035 

0.038 

StdDev 

0.982 

0.086 

0.007 

0.005 

0.014 


Table 1. Registration Accuracy of Various 
Methods 
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Figure 14. Fill Factor Effect 


Figure 14 displays the registration accuracy of Method 4 vs. the FPA fill factor. The noise level 
was set to zero for this study. The fill factor introduces a small error into the frame registration 
process. This point has not been mentioned in the literature.^'^ 

SUMMARY 

A reliable algorithm for subpixel accuracy frame registration is needed to accurately process 
multispectral imagery. Three different extensions of the popular PCM (phase correlation 
method) to subpixel frame registration were evaluated using a common set of satellite images. 
The test images derived from the satellite images include real-world effects such as nonperiodic 
boundary conditions, dead space between pixels, and additive noise. The results are as follows: 

1. Algorithm Comparison: 

• Method 3 (GTF) and its two minor variants (Methods 4 and 5) performed best, 
with registration errors consistently on the order of 0.05 pixels or less. This 
registration accuracy pertains to 256 x 256 images with rms noise levels as high 
as 10%. 
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• Method 2 (RVJ) performed inconsistently. It worked as well as Method 3 for most 
images, but poorly for others. 

• Method 1 (Hoge) is not recommended in its present form. However, the 
mathematical basis is sound, and with minor changes it might be made to 
perform as well as the other methods. 

2. Insignificant Gains from Preprocessing: 

Attempts to further increase the accuracy of registration by preprocessing the images 
(e.g., by taking gradients or by performing histogram equalization of the intensities) led 
to only minor gains in accuracy. 

3. Effect of Focal Plane Array (FPA) Fill Factor 

The fill factor of the FPA sets a limit upon the achievable accuracy. However, this 
effect will generally be only a small part of the error budget for most imaging systems. 
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